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Quantitative Characterization of the Micro structure and Transport Properties of Biopolymer Networks2 

Abstract. 

Biopolymer networks are of fundamental importance to many biological 
processes in normal and tumorous tissues. In this paper, we employ the 
panoply of theoretical and simulation techniques developed for characterizing 
heterogeneous materials to quantify the microstructure and effective diffusive 
transport properties (diffusion coefficient De and mean survival time t) of collagen 
type I networks at various collagen concentrations. In particular, we compute 
the pore-size probability density function P{S) for the networks and present a 
variety of analytical estimates of the effective diffusion coefficient De for finite- 
sized diffusing particles, including the low-density approximation, the Ogston 
approximation, and the Torquato approximation. The Hashin-Strikman upper 
bound on the effective diffusion coefficient Dg and the pore-size lower bound on the 
mean survival time t are used as benchmarks to test our analytical approximations 
and numerical results. Moreover, we generalize the efficient first-passage-timo 
techniques for Brownian-motion simulations in suspensions of spheres to the case 
of fiber networks and compute the associated effective diffusion coefficient De as 
well as the mean survival time t, which is related to nuclear magnetic resonance 
(NMR) relaxation times. Our numerical results for De are in excellent agreement 
with analytical results for simple network microstructures, such as periodic arrays 
of parallel cylinders. Specifically, the Torquato approximation provides the most 
accurate estimates of De for all collagen concentrations among all of the analytical 
approximations wc consider. We formulate a universal curve for r for the networks 
at different collagen concentrations, extending the work of Ycong and Torquato 
[J. Chcm. Phys. 106, 8814 (1997)]. We apply rigorous cross-property relations to 
estimate the effective bulk modulus of collagen networks from a knowledge of the 
effective diffusion coefficient computed here. The use of cross-property relations 
to link other physical properties to the transport properties of collagen networks 
is also discussed. 
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1. Introduction 

Biopolymer networks, such as the cross- hnked bundles (or fibers) of collagen and fibrin 
in the extracellular matrix (ECM), provide mechanical support for cells and serve as 
the media for the transmission of many biomechanical/biochemical cues that regulate 
cell motility, proliferation, differentiation and apoptosis [U [2l [3l H] ■ The diffusion and 
absorption of various macromolecules in biopolymer networks are of crucial importance 
to the regulation and metabolism of normal organs and to the delivery of drugs 
in tumor tissues [SJ [S]- Such biological processes are largely determined by the 
composition and microstructure of the network, especially the complex pore space 
between the fibers [71 H]. Thus, knowledge of the effective transport and mechanical 
properties of biopolymer networks is crucial in order to understand quantitatively the 
aforementioned biological processes. 

The microstructure of biopolymer networks and their associated transport 
properties have been investigated by many researchers. For example, Ogston et al. 
[HI Uni introduced the idea of "influence cylinders" associated with each fiber in the 
system, which enables one to obtain the probability distribution of spherical pores 
with different sizes J, i.e., the pore-size probability density function P{S) (also referred 
to as the pore-size distribution function in literature; see the definition in Sec. 2). For 
polymer networks composed of very long and stiff fibers, Ogston derived an analytical 
expression of -P((5), which depends on the volume fraction and the characteristic 
diameter of the fibers [3]. Other approaches used to ascertain pore-size statistics 
include Fourier analyses of three-dimensional (3D) confocal microscopy images |11| . 
or statistical analysis of nearest points on collagen fibers obtained from confocal- 
microscopic image stacks [12]. Recently, a method based on the direct analysis of 
the entire 3D real-space network geometry from high-contrast confocal microscopy 
data has been developed [131 [T3] . Specifically, Lindstrom et al. [13] represented the 
collagen fibers in the networks as thinned skeletal center lines and the cross-links 
are represented as nodes, which can be thought of as the "graph" representation 
of a biopolymer network. These authors also employed an inverse reconstruction 
method to characterize the microstructure of the collagen networks and investigated 
the mechanical properties of the networks using finite-element analysis. 

The determination of the effective diffusion coefficient D^ for polymer networks 
dates back to the pioneering work of Johansson, Lofroth and coworkers [TS] [TBI HZl 
[T8l [T9] . Johansson et al. experimentally studied the diffusion of small monodisperse 
polyethylene glycols J15j and nonionic micelles |16j in polymer systems and accurately 
measured the long-time- limit self-diffusion coefficient (i.e., D,,) using a tracer technique 
[T7] . By considering the local diffusion of a particle around a single fiber, Johansson 
et al. |18| derived an analytical approximation of D^ which incorporates the 
microstructural information of the pore space. Since their approach was based on 
the key concept of the distribution of "infiuence cylinders" introduced by Ogston, 
the approximation of D,, is henceforth referred to as the Ogston approximation. To 
test the predictive capacity of their theory, Johansson and Lofroth [TO] carried out 
hard-sphere Brownian-motion simulations, in which the diffusing particles were hard 
spheres and the fibers were considered to hinder the diffusion of the particles. Although 
hydrodynamic effects were not taken into account in their simulations, the results were 
shown to be in excellent agreement with experimental data and theoretical predictions 
for a wide range of particle sizes [H] . Recently, the effects of the anisotropy of fiber 
orientations [20] [21] and of the hydrodynamic interactions between the particle and 
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fibers [H] on the effective diffusion coefficient liave also been investigated. We note 
that by mathematical analogy, the problem of macromolecular diffusion in the pore 
space exterior to the collagen fibers is equivalent to the electrical or thermal conduction 
problem in the pore space with perfectly insulating fibers [23] . 

Very recently, it has been suggested that the powerful theoretical and simulation 
techniques developed for characterizing the microstructure and effective properties 
of random heterogeneous materials [23l [24l ESI ES] could be fruitfully employed to 
model complex biological systems, especially malignant tumors and the associated host 
microenvironment [37] . This idea has led to fruitful applications in the understanding 
of the spatial organizations of abnormal cells in brain tumors |28| . 

In this paper, we further explore these techniques from the theory of 
heterogeneous materials by investigating the microstructure and transport properties 
of collagen type I networks (i.e., the most abundant collagen of the human body 
found in tissue and bones, and therefore, called "type I"; see Fig. 1). There exist 
analytical expressions that relate the effective transport and mechanical properties 
of general heterogeneous materials to their microstructure via a variety of n-point 
correlation functions; see Ref. [23] and references therein. This formalism has led to 
the evaluation of effective transport and mechanical properties for a variety of classes of 
microstructures, including dispersions of penetrable [29ll3Qll3l], impenetrable spheres 
[35] [33] [31] , oriented fibers [35] [35] and ellipsoid suspensions (37] |23] , fluid-saturated 
rock j38j . and interpenetrating ceramic-metal composites |39| . 

In the case of collagen-like networks, we calculate for the first time a variety 
of structural descriptors as well as the associated transport properties, such as the 
effective diffusion coefficient D^ and mean survival time r of a Brownian particle 
assuming that the fiber interface is perfectly absorbing (i.e., the average time that a 
Brownian particle spends in the solvent before it gets trapped by the fibers). The 
mean survival time is related to the nuclear magnetic resonance (NMR) relaxation 
times as discussed below. The latter transport property is intimately related to 
the pore statistics [23] [40] . We also employ a variety of approximation schemes for 
the effective diffusion coefficient D^, including the low-density approximation |23j . 
Ogston approximation J18j and the Torquato approximation based on the perturbation 
(phase- property contrast) expansion of D^ j41i 123] . These approximations incorporate 
different levels of microstructural information in terms of various lower-order 
correlation functions that statistically describe the network microstructure. The 
Hashin-Strikman upper bound on the effective diffusion coefficient D^ |42| and the 
pore-size lower bound on the mean survival time r [401 143] are used as benchmarks to 
test our analytical approximations as well as numerical simulations. Specifically, we 
generalize the efficient first-passage-time techniques for Brownian-motion simulations 
in suspensions of spheres [44] [45] [46] [47] [48] to the case of fiber networks and 
compute the associated effective diffusion coefficient De and mean survival time r. Our 
numerical results of I?e are in excellent agreement for the analytical results of simple 
network microstructures such as periodic arrays of parallel cylinders. Moreover, we 
show that the Torquato approximation provides the most accurate estimate of D,, for 
all concentrations among the employed approximation schemes. We also formulate a 
universal curve for r for different networks at different collagen concentrations, i.e., we 
devise a way to scale r in such a way that the scaled data for different collagen networks 
collapse onto a single curve. Rigorous cross-property relations [53] are applied to 
estimate the effective bulk modulus of collagen networks from a knowledge of the 
computed effective diffusion coefficient. 
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The rest of the paper is organized as follows: In Sec. 2, we define the statistical 
descriptors that will be used to characterize the network microstructures. In Sec. 3, 
we provide analytical approximations and rigorous bounds for the effective properties 
De and r. In addition, we discuss the first-passage-time simulation techniques for 
network structures in detail. In Sec. 4, we present the analytical and numerical 
results of the correlation functions and the effective transport properties. In Sec. 5, 
we estimate the effective bulk modulus of collagen networks using the cross-property 
relations from the computed effective diffusion coefficients of the networks. In Sec. 
6, we make concluding remarks, including the use of cross-property relations to link 
other physical properties to the transport properties of collagen networks. 

2. Network Microstructure Characterization 

A collagen network can be considered to be a two-phase heterogeneous material 
composed of a fiber phase and solvent phase (i.e., the pore space), which is exterior to 
the fibers. A important feature of such a network microstructure is that both phases 
percolate across the system, i.e., there is a continuous path between any two point of 
the phase of interest that is entirely in the phase of interest, even when the volume 
fraction of the fiber phase (fraction of space covered by the fibers) is very low. In this 
section, we will introduce various statistical microstructural descriptors for a general 
two-phase material, including the n-point correlation functions S'„ and the pore-size 
probability density function P{S). 

2.1. n-Point Correlation Functions 

Consider a two-phase heterogeneous material in which each phase has volume fraction 
(pi (i = 1, 2), it is customary to introduce the indicator function X^*^(x) defined as 

^•"w={J: III w 

where Vi is the region occupied by phase i and V,; is the region occupied by the other 
phase. The statistical characterization of the spatial variations of the material involves 
the calculation of the standard n-point correlation functions: 

5«(xi,X2, • • • ,x„) = (lW(xi)lW(x2) . • •lW(x„)) , (2) 

where the angular brackets (• • •) denote an ensemble average. The quantity 
Sn {'Xi,X2, . . . ,Xn) also givcs thc probability of finding n points positioned at 
Xi , X2 , . . . , x„ all in phase i. 

For statistically homogeneous materials, the n-point correlation function depends 
not on the absolute positions but on their relative displacements, i.e., 

5'i''(xi,x2,---,x„) = 5^*'(xi2,---,xi„), (3) 

for all n > 1, where x^ = Xj — x^. Thus, there is no preferred origin in the system, 
which in Eq. ([3]) we have chosen to be the point Xi. In particular, the one-point 
correlation function is a constant everywhere, namely, it is equal to the volume fraction 
4>i of phase i, i.e., 

5f' = (l«(x)\=0„ (4) 
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which is the probability that a randomly chosen point in the material belongs to phase 
i. For statistically isotropic materials, the n-point correlation function is invariant 
under rigid-body rotation of the spatial coordinates. For n < d, this implies that 5,\ 
depends only on the distances Xij = |xij| (1 < i < j < n). 

In general, an infinite set of S^ with n ~ 1,2,3... is required to completely 
determine the microstructure and thus, the effective properties of a heterogeneous 
material |23j . Specifically, the effective property of interest can written as an infinite 
series involving integrals of such correlation functions. In practice, a complete 
knowledge of all of the Sn is never available. However, it has been shown that certain 
approximations that can be regarded as resummations of the infinite series expansion 
that incorporate lower-order S'„ (e.g., 52, 6*3 and S^ can provide accurate estimates 
of the effective properties [S^ . We note that since only certain weighted integrals of 
the correlation functions are needed, excellent approximations of effective properties 
can be obtained even without computing all of the Sn explicitly. We will discuss these 
approximations in detail in Sec. 3.1. 

2.2. Pore-Size Probability Density Function 

The pore-size probability density function P{S) first rose to characterize the pore 
space in porous media |43j and was then generalized to characterize any heterogeneous 
material [23]. For a statistically homogeneous and isotropic material, P(d)d6 gives the 
probability that a randomly chosen point in the pore space lies at a distance between 
6 and d + dS from the nearest point on the pore-solid interface. Since it is a probability 
density function with dimensions of inverse length, we have P{S) > for all 5, and it 
normalizes to unity, i.e.. 



P{S)dS = 1. (5) 

/o 

At extreme values of P{6), we have that 

F(0) = s/01, P(oo) = 0, (6) 

where s is the pore-solid interface area per unit volume and i^i is the volume fraction 
of the pore space. Therefore, s/0i is the interface area per unit pore volume. The 
moments of P{5), defined as 



< 5" >= / S''P{S)dS, (7) 

Jo 

provide useful characteristic length scales of the pore space in the material. Certain 

lower-order moments of P{6) also arise in bounds on the mean survival time r |43[ 140) . 

which we will discuss in Sec. 3. 

It is very difficult to obtain analytical expression of P{6) for a general polymer 

network. For networks composed of very long and stiff polymer fibers, Ogston [9] 

derived an expression for P{6), i.e., 

2MS + a) ^^MS+af/a\ (8) 

where (j)2 = 1 — (f>i is the volume fraction of the fibers and a is the fiber radius. For 
(5 = 0, Eq. dS]) gives 

F(0) = s. (9) 
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Comparing Eq. ^ and Eq. ([5]), it is clear that the Ogston expression for P{5) can 
only provide good estimates of it at very low fiber volume fractions, i.e., 02 ^" 
and 01 = 1 — </)2 « 1. In general, the Ogston expression will underestimate P{5) at 
intermediate J, as we will show in Sec. 4. 

Given any network microstructure, the associated P{5) can be numerically 
computed. Specifically, one generates many test points that are randomly distributed 
in the pore space exterior to the fibers and compute the distances from each test point 
to the nearest fiber surface. This amounts to finding the largest test sphere centered 
at the test point that is entirely in the pore space. The resulting distances are binned 
to obtain a probability density function, which is then normalized to give P{5). 

3. Mean Survival Time and Effective Diffusion Coefficient 

S.l. Theoretical Techniques 

3.1.1. Mean Survival Time Consider the steady-state problem of diffusion of 
macromolecules which are absorbed upon contacting the network fibers. This implies 
that the rate of production of the macromolecules per unit volume G is exactly 
compensated by the rate of removal by the traps. Locally, the process is described by 
the following Poisson equation |23 

i:»iAc= -Gin Vi, c = on9V, (10) 

where c is concentration of the macromolecule, Di is diffusion coefficient of the 
macromolecule in the pore space Vi and dV is the pore-fiber interface. The boundary 
condition that c = on dV assumes a perfectly absorbing interface, i.e., a diffusion- 
controlled reaction. This boundary condition can easily be relaxed to take into account 
partially absorbing interfaces |231 140] . 

An important quantity is the mean survival time r associated with a 
macromolecule which is the average time that a diffusing macromolecule spends in 
the pore space before it gets trapped by the fibers. In many medical applications, 
the efficiency of a drug strongly depends on the ability of the drug macromolecules to 
diffuse through the extracellular space mainly composed of collagen without getting 
trapped by the fibers [5]. It is noteowrthy that nuclear magnetic resonance (NMR) 
relaxation in porous media, a widely used technique for biomedical imaging, yields an 
NMR relaxation times from which one can extract the mean survival time we consider 
here [23]. 

The mean survival time r is inversely proportional to the trapping constant 7, 
i.e., 

where 7 is defined via 

G = 7L»i < (12) 

where < c > is the ensemble-averaged concentration field. 

The optimal lower bound on the mean survival time t that incorporates 
information on the pore space in terms of the first moment of the pore-size probability 
density function is given by |401 H5] 

. . ^, (13, 
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where < 5 >= J„ 5P{S)d5. It can be seen from Eq. p^ that tDi can provide an 
estimate of the average pore size of the network. A corresponding upper bound on 
the trapping constant can be obtain by substituting Eq. ([TT]) into Eq. P^ . 

3.1.2. Effective Diffusion Coefficient Consider macromolecules diffusing between 
the fibers that are not absorbed by the fibers. The Brownian motions of the 
macromolecules are hindered by the fibers in the network which results in an effective 
diffusion coefficient D^ smaller than that of the pure solvent in the pore space Z?i. 
By mathematical analogy, the problem of macromolecular diffusion in the pore space 
exterior to the fibers is equivalent to the electrical or thermal conduction problem 
in the pore space with perfectly insulating fibers. Specifically, the local "fiux" J(x) 
of macromolecules is proportional to a local "intensity" E(x) which is the negative 
gradient of the macromolcculc concentration field c(x), i.e., 

J(x) = D(x) • E(x) ^ -D(x) • Vc(x), (14) 

where 

D(x) = I ■^^"^' ^ ^ "^1 (15) 

^ ' \ 0, otherwise 

and I is the unit second-order tensor. Under steady-state conditions with no source 
and sink terms, the conservation of macromolecules requires that J(x) be solenoidal 
[23], i.e., 

V-J(x)=0. (16) 

If D(x) in Eq. (fT4)) is replaced by the local conductivity tensor ct(x), one 

obtains the local governing equations for conduction problems. We would like to 

emphasize that although the diffusion problem and conduction problem are equivalent 

in their mathematical formulations, there is an important distinction between the 

effective diffusion coefficient Dg and the effective conductivity a,,- For the conduction 

problem, although the fiber phase is insulating, its contribution to a,^ is still explicitly 

considered. For example, suppose that one randomly places test particles and tracks 

their Brownian motions to compute a^. There is a fraction of total number of particles 

(f)2, which are initially in the insulating fiber phase and will be trapped there forever. 

Clearly these test particles, which have a zero diffusion coefficient, are taken into 

account in the ensemble average for a^- On the other hand, only the test particles in 

the pore space are considered in order to compute Z^e- Therefore, D^ and Ue for the 

same microstructure are related to one another via the following relation: 

Dp I a,, 1 ap , ^ 

—!- = ^ (17) 

Di 4>iai (l-02)fTi ^ ' 

In the following, we present rigorous bounds and various analytical approxima- 
tions for the effective diffusion coefficient Dp . These results were reported in literature 
for the effective conductivity Up of a general heterogeneous material. Here, we modify 
them according to Eq. ([TT]) to obtain expressions for Dp. 

Hashin-Strikman Upper Bound 

For a two-phase heterogeneous material with an arbitrary but isotropic 
microstructure in which one of the phases (e.g., phase 2) is insulating, the Hashin- 
Strikman (HS) upper bound for the effective diffusion coefficient Dp is given by 

^ = ^. (18) 

Di 2 + 62 ^ ' 
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The HS lower bound in this case is triviaUy zero for all values of 02 [13] • 

Although only volume fractions explicitly appear in the expression, it has been 
shown that HS bound is the optimal bound given the two-point information of the 
isotropic microstructure, i.e., 5*2 (23j . Specifically, it is shown that the HS bounds are 
realizable for a special class of "coated sphere" model microstructirres [23]. However, 
it is clear that such two-point information is far from a complete characterization of 
the network microstructure. Therefore, it can be expected that HS upper bound is 
not tight, as we will show in Sec. 4. 

Low-Density Approximation 

For a heterogeneous material with microstructure composed of well-defined 
inclusions (such as spheres, ellipsoids or cylinders) in a matrix, the effective properties 
of the material can be written as a power series of the volume fraction 4i2 of the 
inclusions [53]. When (/)2 is sufhciently small, i.e., in the low-density limit, truncating 
the power series through the first-order in 02 can provide a reasonable estimate of the 
effective properties of interest. 

For fiber networks, we consider that each fiber is an elongated prolate spheroid in 
the "needle" limit. In such cases, the effective diffusion tensor for dilute suspensions 



(19) 

om. (20) 

For statistically isotropic materials such as suspensions of randomly orientated needles 
in a matrix, the effective diffusion coefficient is the average of the three principal 
components of the tensor De, i.e., 

§^ = l-^'^2 + O(02). (21) 

Physically, Eq. ([^T]) corresponds to the diffusion of Brownian particles in a 
matrix with a single infinitely long fiber, which clearly underestimates Z^e for actual 
biopolymer networks. However, for certain microstructures such as periodic arrays 
of parallel cylinders at low volume fractions, Eq. (j2ip . which we call the low-density 
approximation^ can provide accurate estimates of Z?e, which can be used as benchmarks 
to test our simulation results. 

Ogston Approximation 

An improved approximation for D^ over the aforementioned low-density 
approximation can be obtained if the contributions of multiple fibers are taken into 
account simultaneously. This can be done by using the idea of an "influence cylinder" 
associated with each fiber introduced by Ogston [TU] . Specifically, consider a "coated 
cylinder" with outer radius h and inner radius a (i.e., the radius of the fiber). One 



of orientated needles is 


given by 


De = 


" (i?e)ll 
(i^e)22 

(7^2)33 


where 


Pe)l 


l = (i?e)22= A[l-02 + 
(i^e)33=i^l[l + O(02)] 
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can easily compute the local effective diffusion coefficient Dj^ib) associated with the 
coated cylinder, i.e., 

^^ = ^ = ^ (22) 

D^ l + a2/62 1 + 02(6)' ^ ^ 

where </>2(&) is local volume fraction of the fiber in the coated cylinder. 

Now consider that the global I?e of a network is a weighted average of the local 
Dl(6) for the influence cylinders associated with each fiber, which leads to the relation 

where f{b) is the influence cylinder distribution function [18^. Ogston and coworkers 
assume that the influence cylinders contribute to the global De in the same way they 
contribute to the pore-size probability density function P{5), i.e., 

/•CJO 

P{5)^ f{h)g{h,5)dh, (24) 

J a 

where g(b,d) is the local pore-size distribution associated with a coated cylinder with 
outer radius b given by 

gib,6)^^^^±^H[5-ib-a)] (25) 

and H{x) is the Heaviside step function, equal to unity for a; > and zero otherwise. 
Then /(6) can be obtained by dc- convolution of Eq. (f24|) with a knowledge of P{S), 
either from direct numerical sampling or theoretical considerations. 

Torquato Approximation 

Torquato has derived an expansion of the effective conductivity ag of any 
two-phase heterogeneous materials in terms of the contrast (difference) between 
the conductivities of the individual phases [H]. He showed that the [2,2] Pade 
approximant of the so-called "strong-contrast" expansion, which incorporates up to 
4-point microstructural information involving integrals over S2 , S3 and S'4 can provide 
excellent estimates of cje for a wide range of model microstructures |41j . 

Here, we present the modified 4-point Pade approximation of the effective 
diffusion coefficient D^ for collagen networks, which is henceforth referred to as the 
Torquato approximation, i.e.. 



D._ 1 (i + ig-K2) + (-i + K2-if) 



(26) 



D, 1-02 (l + ig-iC2) + (i + K2 + ig)</'2 ' 

where the parameter C2 is a weighted integral that involves correlation functions Si, 
S2 and 53 of the fiber phase; and the parameter 72 is a weighted integral that involves 
correlation functions Si, S2, S3 and S'4 of the fiber phase. The readers are referred 
to Ref. [41] for detailed discussions of these parameters. Useful rigorous inequalities 
relating C2 and 72 for three-dimensional microstructure are the following |41j : 

- 1 < 72/C2 < 1 - 2C2. (27) 

It is in general nontrivial to compute the 3-point and 4-point correlation functions 
5*3 and S'4, even for relatively simple model microstructures (e.g., dispersions of 
spheres), to obtain exact values of C2 and 72. However, has been shown that 
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simplfied estimates of these parameters based on limited microstructural information 
in conjunction with Eq. ([25]) can lead to excellent approximations for the effective 
properties [HI [53] . Here, we consider the low-density approximation of C2 for a prolate 
spheroid in the needle limit and only keep the leading order term, i.e., C2 = 1/4 |23) . 
The inequalities given in Eq. (j27p then become 

- 1 < 72/C2 < 1/2. (28) 

It has been shown that for dispersions of spheres, 72 = can provide a very accurate 
approximation formula for the effective conductivity of a wide range of sphere volume 
fractions [H]. We will show in Sec. 4 that a proper choice of 72 value that is near 
its lower bound can provide an excellent approximation for the effective diffusion 
coefhcient associated with biopolymer networks. 

3.2. First-Passage-Time Simulation Techniques 

A straightforward way of numerically studying Brownian motion is to simulate the 
exact zig-zag path of a diffusing particle (e.g., see Ref. [H]). However, it is clear that 
this direct approach is not efficient means of obtaining effective diffusive properties, 
since the details of the diffusion paths are averaged out and do not contribute to 
the effective behavior. Moreover, one needs to consider a wide range of step sizes 
associated with each random Brownian jump to extrapolate the results to the case of 
infinitesimal small step size. 

An alternative but much more computationally efficient approach is the first- 
passage-time (FPT) simulation technique introduced by Torquato and coworkers 
[m [45J [461 [ill [ig. The key idea of the FPT approach is not to simulate the details 
of the zig-zag diffusion paths but rather consider the average time that it takes a 
Brownian particle to "jump" directly to a random location on the surface of the 
largest imaginary sphere that is centered at the original position of the particle and 
entirely within the solvent (i.e., the pore phase). The imaginary sphere is referred to 
the "first-passage sphere" (EPS), whose radius is R. It can be shown that the mean 
time Tfps for a Brownian particle, which is initially at the center of the EPS and takes 
a complicated zig-zag path to hit the surface of the EPS is, in three dimensions, given 

by 

Tpps - i?V(2di?i), (29) 

where d = 3 is the space dimension. 

When the particle is very close to the fiber surface, i.e., the distance r from the 
particle centroid to the fiber surface is smaller than a prescribed tolerance A, we 
consider that the particle hits the fiber surface and is reflected back. The EPS in this 
case encloses both the pore phase and fiber phase. Suppose that the EPS centered at 
the Brownian particle centroid possess a radius R, the associated time Tref that the 
Brownian particle takes to hit the fiber surface, be reflected back and hit the EPS can 
be estimated by 



, , R^ V1+V2 



m=0 



,(30) 



where Q < r < R, Vi and V2 are the volume of the pore phase and the volume of the 
fiber phase enclosed in the EPS, respectively, and 

_ {-ir+\2m)\ 3(4m + 3) 

2™+i 22™+i(m!)2 (2m-l)(m-H2)(m-^l)' ^ ' 
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Equation (PH)) was first derived by Torquato and eoworkers [ISJ HH 27] , and it has been 
shown to provide excellent approximation of the exact t^ef for any local geometry 
when r -^ R and R is smaller than the diameter of the fiber. The readers are referred 
to Ref. [35] and the references therein for additional details. 

To compute Tref using Eq. ((30)) . a key step is to evaluate the intersection volume 
between a sphere and a cylinder (i.e., V2), the details of which are given in the 
Appendix. In the rare case that the particle is close to a cross link (junction of several 
fibers), V2 is the volume of the fiber junction enclosed in the FPS, which is computed 
by Monte Carlo sampling [17]. For example, one randomly places test points in the 
FPS and computes the fraction of times that the point falls into the vicinity of the 
fiber junction. 

To obtain D^^, one considers an ensemble of Brownian trajectories in the pore 
space. When a diffusing particle is sufficiently far away from the fiber surface, one 
constructs the largest FPS of radius R around the diffusing particle which just touches 
the fiber surface. The particle then jumps in one step to a random point on the surface 
of the FPS and the process is repeated, each time keeping track oi Rf, until the particle 
is within a prescribed very small distance A to the fiber surface (see Fig. 3). At this 
point in time, the particle is considered to hit the fiber and then is refiected back. 
Thus, one keeps track of the radius Rj of the FPS that encloses both the fiber phase 
and the pore phase and computes the associated time t-i^-ep{Rj). The expression for 
effective diffusion coefficient D^ is then given by 

where Tref {R) is given by Eq. pop and < . > denotes the ensemble average over many 
Brownian particles. In our simulations, we use N — 5 000 Brownian particles. 

The mean survival time r can be obtained in a similar way f44[ . Specifically, one 
constructs the first-passage-time path composed of many jumps to the surface of FPS 
associated with a Brownian particle and keeps track of Ri for each FPS. When the 
particle is within A to the fiber surface, it is considered trapped by the fiber. Thus, 
the mean survival time can be computed via 

T^ly^R^/DiV (33) 

where < . > denotes ensemble average over many Brownian particle. In our 
simulations, we use N = 5 000 Brownian particles. 

We note that in the aforementioned first-passage-time simulation technique, we 
consider that the Brownian particle is "point" particle with zero diameter. For finite- 
sized particles diameter dp, it has been shown that one can still consider "point" 
particles in a network microstructure with the diameter of the fibers dp dilated by dp 



4. Results 

Our data are "graph representations" of collagen type I networks [14] with final 
collagen concentrations of 1.0 mg/ml, 2.0 mg/ml and 4.0 mg/ml. The fibers roughly 
possess a circular cross section of diameter dp = 1.0 x 10~^ m [T3]. The average fiber 
lengths ip for the networks with the three collagen concentrations are respectively 
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1.96 X 10~^ni, 1.81 X 10~^m, and 1.28 x lO^^m. The corresponding volume fraetions 
of the fibers are respectively 1.7 x 10~^, 2.4 x lO"'', and 5.2 x lO^''. The tolerance 
A described in Sec. 3 is chosen to be A = 5 x 10"^^^ = 5.0 x 10~^° m. Since we 
do not consider hydrodynamic effects in our simulations, we only consider particles 
with diameter dp comparable to the fiber diameter dp, i.e., dp < dp. For large dp, 
it has been shown that the hydrodynamic effects on Brownian motion are significant 
[22| . The results reported below are ensemble averages of three independent network 
configurations at each collagen concentration. 

4.I. Pore-Size Probability Density Function 

The pore-size probability density functions P{6) for the collagen network at three 
concentrations are numerically computed as described in Sec. 2. The obtained P{5) 
are shown in Fig. 4 and compared to the corresponding Ogston expressions [Eq. ([5])] 
at the same fiber volume fractions. 

It can be clearly seen that the Ogston expression oi P{6) overestimates the number 
of intermediate pores and underestimates the number of large pores in the system. 
This is because in the derivation of Eq. ([5]), it is assumed that the network is composed 
of fibers with very long persistence length. For the collagen networks we study, the 
average fiber lengths £ are less than twice of the corresponding averaged pore size 
< S > [defined in Eq. ([7])], which are respectively 1.22 x 10~^m, 0.998 x 10~^m, 
and 0.684 x 10~^m for collagen concentrations 1.0 mg/ml, 2.0 mg/ml and 4.0 mg/ml. 
Therefore, the long-fibcr-length assumption for the Ogston expression is not true here. 

The P{S) data will be employed to compute the lower bound on the mean survival 
time r [Eq. ([T3| ] and to compute the Ogston approximation for the effective diffusion 
coefficient D^ [Eq. ([25)1 ] in the following sections. 

4-2. Mean Survival Time 

The mean survival time r is computed using the first-passage-time technique described 
in Sec. 2. Figure 5 shows the scaled dimensionless mean survival time tDi/£'^ for 
the collagen networks with three different concentrations as a function of Brownian 
particle diameter. It can be seen that as the collagen concentration increases, larger 
particles are more easily get trapped by the fibers. This fact is of great importance in 
cancer chemotherapy, which we will discuss in Sec. 6. 

As indicated in Sec. 3, the diffusion of finite-sized particles in the original network 
is equivalent to the diffusion of point particles in a properly dilated network, which 
possesses a higher fiber volume fraction. Figure 6 shows the scaled mean survival time 
tDi/IJ^ for the collagen networks with three different concentrations as a function 
of the particle diameter. The pore-size lower bounds are also shown. It can be 
seen that although the bounds are not sharp, they do not deviate very much from 
the actual mean survival times. We note that these bounds only incorporate partial 
information about the pore-size probability density function P{S), namely, the first 
moment < S > oi P{S). Therefore, one would expect that incorporating the full 
information content of P{6) would lead to good predictions of the effective diffusive 
properties considered here. Indeed, we will show in the following section that the a 
generalization of the Ogston approximation that employs the complete microstructural 
information contained in P{S) provides a very good estimate of D^- 

In Ref. |50] . Torquato and Yeong found a universal curve for a scaled mean 
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survival time r for a wide range of microstructures with different porosities, including 
various random and ordered distributions of spheres and certain continuous models. 
Specifically, the universal curve has the following form: 

— ~aix + a2x'^, (34) 

To 

where oi and 02 are constants and 

To = r, , 2 , a;= -— , (35) 

and s is the specific surface, i.e., the solid-pore interface area per unit volume. For 
the class of microstructures they studied, Torquato and Yeong found that ai = 8/5 
and 02 = 8/7. 

For the collagen networks studied here, we find that Eq. p4p also holds (see 
Fig. 7). However, the constants are different from those obtained by Torquato and 
Yeong, i.e, we find that ai = 0.121 and 02 = 1.88 for the networks with different 
concentrations. A possible reason for the difference in the constants is that collagen 
networks do not belong to the same class of microstructures studied in Ref. [50] , which, 
for example, do not contain filamentary-like structures, as in the case of collagen fibers. 
This implies that there could exist a more general scaling curve for the mean survival 
time that incorporates both the networks and the microstructures studied in Ref. [SO] . 
Nonetheless, our results enable one to efficiently estimate the properties of collagen 
networks. In particular, given any of the three quantities among the four quantities 
T, 4>i, s and < 5 >, the remaining one can be estimated employing Eq. p4p . 

4.-3. Effective Diffusion Coefficient 

The effective diffusion coefficient Dp for various network microstructures are computed 
using both the theoretical techniques described in Sec. 3.1 and the first-passage-time 
technique described in Sec. 3.2. Figure 8 shows D^ for the fiber networks with different 
collagen concentrations as a function of the Brownian particle diameter. Similar to the 
case of the mean survival time, as the collagen concentration (p increases, it becomes 
more and more difficult for larger particles to diffuse in the collagen. 

Figure 9 shows D^ as a function of the fiber volume fraction. In addition to 
the results for the collagen networks, we also show D^, for a model microstructure 
composed of parallel cylinders arranged on a square lattice. The Hashin-Strikman 
(HS) upper bound and various approximations of D,. discussed in Sec. 3.1 are also 
shown in Figure 9. As we indicated earlier, since the HS bound only incorporates 
the limited two-point information 6*2, it can not provide a good estimate of D^. 
This is also evident from the fact that the HS bound is realized by certain class of 
"coated sphere" model microstructures, which are clearly topologically distinct from 
the network microstructures because one of the phases is topologically disconnected. 
It can be seen that the low-density approximation also underestimates D^ for the 
collagen networks. This is because it only considers the effect of a single long fiber 
to the diffusing particles. In the actual networks, the average fiber length is less than 
twice the average pore size as we indicated earlier. Moreover, the cross-links also 
significantly hinder the diffusion of the particles. However, for the parallel-cylinder 
model at low volume fractions, the low-density approximation should provide accurate 
estimates, since in such cases, the diffusion of the particles is only hindered by well 
separated single cylinders. Indeed, we find that the approximation agrees very well 
with our simulation data, which also verifies the accuracy of our simulations. 



Quantitative Characterization of the Micro structure and Transport Properties of Biopolymer Networkslb 

The Ogston approximation that incorporates the porc-sizc information of the 
networks provides a much better estimate of Dp compared to the HS bound and the 
dilute approximation. However, it still slightly underestimates D^. for large particles 
in networks at high collagen concentrations. This is because the "influence cylinders" 
(see Sec. 3.1) are associated with individual long fibers and the effects of finite-fiber 
length and the cross-links are still not fully incorporated. On the other hand, one can 
see that the Torquato approximation agrees extremely well with the simulation data 
for all volume fractions that we considered. Specifically, in employing Eq. (|26|) we have 
chosen the 4-point parameter value such that 72/C2 = —0.925. Note that this value is 
very close to the lower bound value -1, which is not very surprising the networks can 
be considered as a kind of "limit" microstructure. This does not mean the actual value 
of 72/C2 if computed with a full knowledge of the associated 3-point function S^ and 4- 
point function S'4, since the value of C2 we used is also an approximation. The success 
of the Torquato approximation is due to the fact that higher-order microstructural 
information that possibly reflects the effects of the cross-links is already taken into 
account by the 4-point parameter 72 in the expression Eq. (|26|) . 



5. Estimating Elastic Properties of Collagen Network Using 
Cross-Property Relations 

Since effective properties of heterogeneous materials reflect certain microstructural 
information about the material, it is possible to extract rigorously information about 
one physical property given an accurate determination of a different effective property 
obtained either experimentally or theoretically. Such interrelationships are called 
cross-property relations |23[ 1561 I57[ 1581 1591 160] . Rigorous cross-property relations 
become especially useful if one property is more easily measured than another property. 
In this section, we estimate the effective bulk modulus K^ [13 of fluid-saturated 
collagen networks using the effective diffusion coefficient De computed here using the 
cross-property relations. In particular, Gibiansky and Torquato [Ml [57l [58] derived a 
nontrivial rigorous cross-property upper bound K^ on the effective bulk modulus Ke 
of a fluid-saturated porous material with an insulating solid phase given the effective 
conductivity (equivalent to the effective diffusion coefficient) of the material, i.e., 

K. < K - ^le - „[3,^^f_3„_202G2]' ^^^^ 

where 

a = 02i^i + 0ii^2 + 4G2/3, 
Ki, = 0iXi + 02i^2 - 0i02(A'i - K2)ya, 

and 01 and Ki are respectively the volume fraction and bulk modulus of the fluid 
phase; 02, K2, G2 are respectively the volume fraction, bulk modulus and shear 
modulus of the solid phase. Moreover, the "formation factor" F in Eq. ([36]) is given 

by 

1 ^1 

01 De 

where Di is the diffusion coefficient of the fluid phase and D^ is the effective diffusion 
coefficient of the porous material. 

For the collagen networks considered here, the solid phase corresponds to the 
collagen fibers. The bulk modulus Kg. of fluid-saturated collagen networks at fiber 



2 /. (37) 
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volume fraction 02 — 0.005 has been measured experimentally [ST], i.e., K^ w 2500 
Pa. The shear modulus of "dry" collagen networks (i.e.., fiber networks without fluid) 
have been computed numerically by Lindstrom et al. [TJ, i.e., G2 = 24 Pa. We 
use K2 = 10 Pa for the "dry" network and A'l — 2 GPa for the fluid. Using the 
Torquato approximation for the effective diffusion coefficient D^, the upper bound 
value obtained from inequality (|36|) is computed, i.e., K^/ = 3530 Pa, which provides 
a surprisingly good estimate of Ke, given the fact that K^^ is a rigorous upper bound 



6. Conclusions and Discussion 

In this paper, we have quantitatively characterized the microstructure, the mean 
survival time r, and the effective diffusion coefficient D^ of collagen type I networks by 
applying theoretical and computational techniques from the theory of heterogeneous 
materials. Specifically, we have computed the pore-size probability density function 
P{S) for the networks. We have also employed a variety of theoretical approximation 
schemes for the effective conductivity of a two-phase material to estimate the effective 
diffusion coefficient De for the networks. Such estimates include the low-density 
approximation, the Ogston approximation, and the Torquato approximation, all of 
which incorporate different levels of microstructural information about the networks. 
The Hashin-Strikman upper bound on D,, and the pore-size lower bound on r are used 
as benchmarks to test our results. Moreover, we have generalize the efficient first- 
passage-time techniques for Brownian-motion simulations in suspensions of spheres to 
the case of network microstructures and compute the associated Df, and r. We have 
found a universal curve for r for the networks at different collagen concentrations 
and have shown that the Torquato approximation which takes into account higher- 
order microstructural information can provide the most accurate estimate of Dg 
for all collagen concentrations among the employed approximation schemes. Our 
work also demonstrates that employing the rich family of theoretical and simulation 
techniques developed in materials sciences to characterize biological systems (e.g., 
the heterogeneous host microenvironment of tumors) suggested in Ref. [27] is a very 
promising approach worthy further exploration. 

We have found that as the collagen concentration increases, the diffusion of large 
particles in the collagen networks, and thus the extracellular matrix (ECM) becomes 
more and more difficult and it is easier for the diffusing particles to be trapped by the 
fibers. This is a major problem associated with any cancer chemotherapy, since drug 
macromolecules would get trapped by collagen fibers without successfully diffusing to 
the target site. It is known that a growing malignant tumor constantly modifies the 
chemical composition of the collagen networks composing its ECM [51] . In addition, 
since a pressure is built up as the tumor grows, the surrounding ECM is pushed 
and compressed, leading to a higher collagen concentration in tumor ECM than in 
normal tissues [551 [531 [SI] . Therefore, it can be expected that the diffusion of drugs to 
the tumors is really difficult. More efficient chemotherapies trying to overcome these 
difficulties are being developed [55] . 

We also applied a rigorous cross-property upper bound to to estimate the effective 
bulk modulus K^ of collagen networks from a knowledge of the effective diffusion 
coefficient D^ computed here. The estimated value of Kg agrees well with existing 
experimental data, given the fact that it is a rigorous upper bound. 

In future work, we intend to generalize our simulation techniques and theoretical 
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approaches to investigate transport properties of tissues with both collagen networks 
and various types of cells. Specifically, we will focus on the effects of the cell shape and 
the plasma membrane on the diffusion of maeromolecules. In addition, we will model 
the mechanical behavior of tissues using the well-developed methods for heterogeneous 
materials. Progress in these studies should deepen our understanding of the effects of 
the host microenvironment on tumor growth and would lead to better cancer treatment 
strategies. 

In addition, we will apply cross-property relations to estimate other physical 
properties of collagen networks from a knowledge of the effective diffusive transport 
properties computed in this paper. In particular, given the latter, we will bound the 
fluid permeability [59l [60] for the collagen networks studied here. 
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Appendix: Intersection volume of a sphere writh a cylinder 

Consider a sphere with radius Rs centered at the surface of a cylinder with radius Re, 
the intersection volume Vi of the sphere and the cylinder for the case Re > Rs is given 

by [ng 

Vi = ^ttRI+^^ [K{k){A - B){3B - 2A) + E{k)A{2A - 45)] , (39) 

where K{k) and E{k) are elliptic integrals of the first and second kind, respectively, 
i.e., 

K{k)= f ^=^=== E{k)^ fdzJ^^^, (40) 
Jo sj [l - z^){l - k^ z^) Jo V 1-2^ 

and 

A = ARl, B = Rl, fc2 ^ 5/^ (42) 
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Abbreviations list 

• ECM: extraccUuar matrix 

• HS: Hashin-Striknian 

• NMR: nuclear magnetic resonance 

• FPS: first-passage sphere 
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Figure Legends 

• Figure 1: Collagen networks, (a) Confocal microscope image of collagen-I 
at a final collagen concentration 2.0mg/ml. The linear size of the image is 
approximately 150 /im. Image courtesy of S. B. Lindstrom. (b) Three-dimensional 
"graph" representation of the collagen network studied here. The linear size of 
the box is approximately 100 fim. 

• Figure 2: An illustration of sampling P{S) of the collagen network in two 
dimensions. Shown are the thin fibers that can possibly intersect as well as a 
test point (red) and its associated largest sphere (blue) entirely within in the 
pore space exterior to the fibers. 

• Figure 3: An illustration of the first-passage-time simulation technique in two 
dimensions. Shown are the thin fibers that can possibly intersect and several first 
passage spheres. Starting from an initial position, a diffusing particle jumps to a 
random location on the surface of its associated first-passage sphere. The process 
is repeated until the particle is very close to the fiber surface. The the particle 
will hits the fiber and be reflected back to the pore space. 

• Figure 4: The pore-size probability density function P{6) for collagen networks 
at different collagen concentrations. 

• Figure 5: The scaled dimensionless mean survival time rDi/Pp for the collagen 
networks at different collagen concentrations as a function of the diffusing particle 
diameter. 

• Figure 6: The scaled dimensionless mean survival time rDi/Pp for the collagen 
networks at different collagen concentrations as a function of the diffusing particle 
diameter and the associated pore-size lower bound. 

• Figure 7: Universal curve for the scaled mean survival time t /tq versus < 5 >^ 
/{tqDi) for the collagen networks at different collagen concentrations. It is seen 
that the scaled mean survival time for different collagen networks collapse onto a 
single curve. 

• Figure 8: The dimensionless effective diffusion coefficient De/Di for the collagen 
networks at different collagen concentrations as a function of the diffusing particle 
diameter. 

• Figure 9: The dimensionless effective diffusion coefficient De/Di for the 
collagen networks at different collagen concentrations as a function of the fiber 
volume fraction. The Hashin-Strikman upper bound and various analytical 
approximations as discussed in Sec. 3.1 are shown and compared to the simulation 
data. 
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Glossary 

• Collagen: a group of naturally occurring and the most abundant proteins 
(biopolymcrs) found in animals, especially in the flesh and connective tissues 
of mammals. 

• Heterogeneous material: a material composed different materials (e.g., a 
composite) or the same material in different state (e.g., a polycrystal). The fluid- 
saturated collagen networks studied here are special heterogeneous materials. 

• Diffusion coefficient: a proportionality constant between the molar flux due to 
molecular diffusion and the gradient in the concentration of the species or the 
driving force for diffusion. 

• Mean survival time: the average time that a diffusing molecule spends in the 
fluid phase before it gets trapped at the interface of the collagen fibers assuming 
a perfectly absorbing interface. 

• Bulk modulus: a measure of a material's resistance to uniform compression, 
defined as the ratio of the infinitesimal pressure increase to the resulting relative 
decrease of the material's volume. 

• Shear modulus: a measure of a material's resistance to shape deformation shear 
stress, defined as the ratio of shear stress to the shear strain. 

• Correlation function: a statistical descriptor of the microstructure of a 
heterogeneous material, quantifying the spatial correlations at different points 
in the microstructure. 

• First-passage-time: the average time that it takes for a diffusing particle to 
"jump" directly to a random location on the surface of an imaginary sphere that 
is centered at the original position of the particle and entirely within the solvent 
region. 

• Cross-property relation: an interrelationship that enables one to extract rigorously 
information about one physical property of a heterogeneous material given an 
accurate determination of a different property. 
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Figure 1. Jiao and Torquato 




Figure 2. Jiao and Torquato 
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Figure 3. Jiao and Torquato 




Figure 4. Jiao and Torquato 
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Figure 5. Jiao and Torquato 
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Figure 6. Jiao and Torquato 
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Figure 7. Jiao and Torquato 
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Figure 8. Jiao and Torquato 
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Figure 9. Jiao and Torquato 



